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We develop a uniform semiclassical trace formula for the density of states of a three-dimensional 
isotropic harmonic oscillator (HO), perturbed by a term ier 4 . This term breaks the U(3) symmetry 
of the HO, resulting in a spherical system with SO (3) symmetry. We first treat the anharmonic 
term for small e in semiclassical perturbation theory by integration of the action of the perturbed 
, periodic HO orbit families over the manifold CP 2 which is covered by the parameters describing 

their four-fold degeneracy. Then we obtain an analytical uniform trace formula for arbitrary e which 
in the limit of strong perturbations (or high energy) asymptotically goes over into the correct trace 
£N| , formula of the full anharmonic system with SO (3) symmetry, and in the limit e (or energy) — > 

restores the HO trace formula with U(3) symmetry. We demonstrate that the gross-shell structure 
O ' of this anharmonically perturbed system is dominated by the two-fold degenerate diameter and 

circular orbits, and not by the orbits with the largest classical degeneracy, which are the three-fold 
degenerate tori with rational ratios u) r : u) v — N : M of radial and angular frequencies. The same 
holds also for the limit of a purely quartic spherical potential V(r) oc r 4 . 
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I. INTRODUCTION 



The semiclassical quantisation of non-integrable systems using properties of their periodic orbits was triggered by 
M. Gutzwiller [| and extended by several groups [H Qi 0| . It allows one to express the oscillating part Sg(E) of 
^vq ■ the quantum-mechanical density of states, given exactly in terms of a quantum spectrum E n and separated into two 
' terms by 

o : 

g(E) = J2S(E-E n )=g(E)+Sg(E), (1) 

in , 

. through a semiclassical trace formula of the form 

m ■ 

Sg sc {E) ~ A po {E) cos [S po (E)/h - a po n/2] . (2) 



The sum is over all periodic orbits (po) of the classical system, S po (E) = § p • dq are their action integrals, the 
amplitudes A po (E) depend on their stabilities and degeneracies, and a po are some phases called Maslov indices. The 
average part g (E) of the density of states, which by definition varies smoothly with energy, is obtained by the extended 
Thomas-Fermi (ETF) model (see, e.g., 6], Chapter 4). Then, the sum gETF^E) + Sg sc (E) usually turns out to be a 
good approximation to the exact quantity QJ , although the sum over the po in J5J is only an asymptotic one, correct 
to leading order in 1/H, and in chaotic systems is hampered by convergence problems Q. For integrable systems, 
Berry and Tabor Q and later Creagh et al. p| showed how a trace formula of the form J5J can quite generally be 
obtained from the Einstein-Brillouin-Kellcr (EBK) quantisation 9]. This method will be used in the present paper 
for a spherical system. 

The semiclassical theory, often referred to as periodic orbit theory (POT), has been very successful not only in 
(partial) quantisation of a given Hamiltonian, but also in interpreting many experimentally observable quantum 
oscillations in finite fermion systems - so-called "shell effects" - in terms of classical mechanics. Examples are atomic 
nuclei (lC|. metallic clusters |Ii1IFaI13| . semiconductor quantum dots |l3L and metallic nanowires 15]. An overview 
over many aspects of the POT and illustrative applications are given in 

One problem that remains with all the trace formulae developed in the work quoted above is that the amplitudes A po 
diverge in situations where a continuous symmetry is broken or restored under the variation of a system parameter 
(which may also be the energy E), or where bifurcations of periodic orbits occur. In such situations one has to 
go beyond the stationary phase approximation which is underlying the semiclassical approach. This has, besides 
El HI EH IM 13 > been developed most systematically by Ozorio de Almeida and Hannay ^(| for both situations, leading 
to local uniform approximations with finite amplitudes A po - So-called global uniform approximations, which yield 



2 



finite amplitudes at symmetry-breaking and bifurcation points and far from them go over into the standard (extended) 
Gutzwiller trace formula, were developed for the breaking of U(l) s ymm etry in for some cases of U(2) and SO(3) 
symmetry breaking in [Tsj , and for various types of bifurcations in . (Details and further references may be found 
in 0, Chapter 6.3.) 

In this paper we investigate an anharmonically perturbed three-dimensional spherical harmonic oscillator (HO) 
Hamiltonian (for a particle with mass m = 1): 

H(r,p) = l -p 2 + \^r 2 + \er 4 = H (r,p) +eSH(r) , (3) 

where p = |p| and r = |r| are the absolute values of the three-dimensional momentum and radius vectors, e > is 
the strength of the perturbation which first is assumed to be small, but later may assume arbitrary positive values. 
The unperturbed HO system has U(3) symmetry [2(| which is broken by the anharmonic term, leading to the SO (3) 
symmetry of the perturbed spherical system. 

The Hamiltonian suggested itself in a recent study |2l| of harmonically trapped fermionic atoms with a short- 
range repulsive two-body interaction, treated self-consistently in the Hartree-Fock approximation. As a result, very 
pronounced shell effects in the total energy E to t of the interacting s yste m as a function of the number N of atoms 
were found, which remind about the so-called super-shells predicted \H\ and observed 0] m metallic clusters. In a 
first attempt to interpret these shell effects semiclassically . 2J| , we parameterised the self-consistent mean field of the 
interacting system by the Hamiltonian @ and applied the semiclassical perturbation theory of Creagh |22j to explain 
qualitatively the shell structure of the HF results. 

In the present paper we describe some of the mathematical details of the semiclassical perturbation theory for 
the symmetry breaking U(3) — ► SO(3) and develop a uniform trace formula valid for arbitrary strengths e in the 
Hamiltonian ©. In Section [H] we present the perturbative trace formula for the density of states which is valid in the 
limit of small e and has been presented shortly in [2l|. It already puts the dominance of the shortest periodic orbits, 
namely the circles and diameters, into evidence. Although their individual contributions to the perturbative trace 
formula diverge in the limit e — > 0, their sum restores to the unperturbed HO trace formula with U(3) symmetry in 
this limit, which at the same time is the limit of zero energy. 

In Section IIHI we develop the uniform trace formula that includes the contributions of the diameter and circle 
orbits for arbitrary values of e. For this purpose, we start from the EBK quantisation of an arbitrary system with 
spherical symmetry and apply the Poisson summation formula. We obtain a one-dimensional trace integral for the 
semiclassical density of states, whose end-point contributions correspond to the diameter and circle orbits, valid for 
an arbitrary spherical one-body potential. For the Hamiltonian ©, the gross-shell structure of the density of states 
is at low energies always dominated by the families of diameter and circle orbits, although these orbit families only 
have a two-fold degeneracy. At sufficiently high energies and repetition numbers, "rational tori" with frequency ratios 
u> r : uo^ — N : M > 2 with N > 7 bifurcate from the circle orbits with repetition numbers M > 3, as discussed in 
Section Till Fl These rational tori have the highest (three- fold) degeneracy possible in a three-dimensional spherical 
system, but due to their length they only affect the finer quantum structures of the density of states. This situation 
is completely different from that of a spherical billiard where the contribution of the diameter orbits to the density 
of states is practically negligible. 

In the limit e — * oo, becomes a purely quartic oscillator which in many respects is easier to handle. In particular, 
in this limit the system acquires the "scaling property" that its classical dynamics becomes independent of the energy 
which can be absorbed by a rescaling of coordinates and time. This limit will be discussed separately in Section 
IIII Gl One interesting classical aspect is that no bifurcations of the circle orbits occur. The same rational tori, which 
bifurcate from the circle orbit in the perturbed HO system, exist here at all energies. But, again, they affect only 
the finer quantum structures of the density of states, while its gross-shell structure is largely dominated by the circles 
and diameters. 

In the appendices 0- 10 we collect some mathematical details about the integration over the manifolds CP 2 and S 5 
and some explicit analytical expressions for action integrals and periods in terms of elliptic integrals. Finally, in the 
appendix [D] we re-derive from our general trace integral the known trace formulae for two popular spherical systems: 
the spherical billiard and the Coulomb potential. 
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II. PERTURB ATIVE TREATMENT OF THE ANHARMONICITY 



We first write down the exact density of states of the unperturbed three-dimensional harmonic oscillator (HO), 
given by the Hamiltonian 

H (r,p) = ±p 2 + ^ 2 r 2 . (4) 
Its quantum-mechanical density of states can be written in the form 

oo 

g (E) = d nS(E - E n ) = g$ TF (E) + Sg (E) , (5) 

using the spectrum E n — hu(n + 3/2) with the degeneracy factor d n = (n+ l)(n + 2)/2. Here g FTF {E) is the smooth 
part given by the ETF model 



9E TF {E) ' 



E 2 -\{h,) 2 



2{hwf 

and Sgo(E) is the oscillating part given by the following exact trace formula 



(6) 



g (E) = 2g FTF (E)J2(~l) k cos ( 2irk^-) = 2 g FTF (E) f>l) fc co S [kS (E)/h] . (7) 



k=l v 7 k=l 

It can be understood as a sum over all periodic orbits of the system. Hereby k represents the repetition number of 
the primitive classical periodic orbit which has the action So(E) — 2ttE/lo. Keeping only the leading TF term of 
g FTF (E) in JfjJ), the density of states can also be written in the form 

9o(E) = ^3 Ke£(-l) fc e ^E/H^ + ^ (8) 

since the imaginary parts cancel upon summation. Note that the smooth TF part in (jfJJ comes from the contribution 
with k = in (0). 

Next we follow Creagh |22j | by writing the perturbed trace formula in the form 

9 P ert(E) = ^3 $teJ2(-l) k M(ka/h)e^ kE /^, (9) 

where A4(ka/h) is a modulation factor that takes into account the lowest-order perturbation of the action of the HO 
orbits. Here a = a(e, E) is a small action that, quite generally, depends on e and the energy E. The factor in front of 
the sum in © takes into account only the leading term of the smooth part g FTF (E); this is consistent with the fact 
that the perturbation theory only deals with the terms of leading order in hT 1 of the semiclassical trace formula. As 
we will see in Section UlI El the contributions of terms of relative order H 2 in g FTF (E) are, indeed, negligible. 

As described in detail in |22j, one obtains M(ka/h) from an integration of a perturbative phase function over the 
manifold that describes the classical degeneracy of the HO orbits. The unperturbed HO has U(3) symmetry, so that 
M.(kajK) is formally obtained by the average 

M(ka/h) = (e lfcAS(o)/? W(3), (10) 

where o is an element of the group U(3) characterising a member of the unperturbed HO orbit family, and AS(o) is 
the lowest-order primitive action shift brought about by the perturbation eSH(r) in J3J). In the present case AS*(o) 
is nonzero already in first-order perturbation theory and therefore a is proportional to e. 

We need, however, not integrate over the full nine-dimensional space of group elements of U(3). It is sufficient to 
consider only the subset of symmetry operations which transform the orbits within a given degenerate family into each 
other (without changing their actions). The dimension of this subset is the degree of degeneracy / of that family. For 
the three-dimensional spherical HO we have / = 4. This can most easily be seen by the following argument which also 
will allow us to find a suitable parametrisation for the four-fold integration. The full phase space is six-dimensional; 
due to energy conservation it is reduced to a five-dimensional energy shell which has the topology of a five-sphere S 5 . 
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Of the five parameters that specify a point on S 5 , one can be chosen as a trivial time shift along the periodic orbits 
corresponding to a simple phase factor e~ luJt . This parameter forms a subgroup U(l), so that its elimination restricts 
us to the four-dimensional manifold (see Appendix 1X1 for its mathematical definition) 

S 5 /U(l) = CP 2 , (11) 

which is neither a four-torus nor a four-sphere [2^ . (Note that for the two-dimensional HO one is led to the manifold 
CP 1 which happens to be homomorphous to the two-sphere S 2 , see |2 2j.) A suitable parametrisation of CP 2 is given 
in Appendix lAl in terms of four angles (i?, ip, i% ^3) whose meaning will become clear in a moment. The modulation 
factor therefore becomes 

M{ka/h) = \l dVL CV 2 e ikAS / H = A f 2 CO s psin pdip f ' sin 3 z9 costf d<d [ dv 2 f dv 3 e ikAS ^'" 2 ' U3 ^ h . (12) 
n J 77 Jo Jo Jo Jo 

In the zero-perturbation limit e — ► 0, where the action shift AS* and hence also a is zero, the modulation factor 
becomes unity, as it should for © to approach JSJ). 

Next we have to determine the action shift AS of a primitive periodic orbit of the HO, caused by the perturbation 
er 4 /4. The harmonic solutions of the classical equations of motion of the unperturbed system are well known: 

Xi = \/2Ei/uj 2 cos(ujt + Vi) , i = l,2,3 {x u x 2 , x 3 ) = (x, y, z). (13) 

Hereby Ei are the three conserved energies in the three dimensions. Depending on the values of Ei and the phases 
v%, i|13|) describes circles, ellipses or librations through the origin; we will in the following call the latter the "diameter 
orbits" . We now re-parameterise the Xi in the following way: 

x(t) = Rn\ cos(wi) , y (t) = Rn 2 cos(w£ + v 2 ) , z(t) — Rn 3 cos(ivt + u s ) , (14) 

where v 2l 1/3 G [0, 27r), and R is given by 

R = V2E/lj , E = E 1 +E 2 +E 3 . (15) 

Since the (ni, n 2l n 3 ) must lie on the sphere S* 2 , due to the conservation of energy, we can also write them as: 

ni=cos$, n 2 — sin i? cos (p , n 3 — sin & sin ip . (16) 

Actually, since the ni must be restricted to positive definite values in order to cover once all classically allowed values 
of the Xi(t), we only have to integrate over the first octant of S 2 , so that tp, $ S [0, w/2]. The four angles $, ip, v 2 , v 3 , 
with their ranges of definition, are precisely the parameters describing the manifold CP 2 as explained in Appendix lAl 
and the correct integration measure is that used in p2|l. We might have kept the phase angle v\ of x{t) in l|14|l . too; 
integration over it is tantamount to a trivial integration over time along the orbits. The resulting five-dimensional 
integral becomes equivalent to that over the S 5 sphere discussed in Appendix iBl 

According to Creagh the first-order action shift brought about by a perturbation eSH is given by 

A 1( S* = -e I SH(t)dt, (17) 

J po 

where po stands for a particular member of the unperturbed periodic orbit family, and the perturbating Hamiltonian 
8H(t) = 5H{xi{t)) has to be evaluated along the periodic orbit. Inserting (fT4^) - (tTB|l into the perturbation 5H = 
r 4 /4 = (x 2 + y 2 + z 2 ) 2 / A leads to elementary integrals over powers of the trigonometric functions. The result is 

AxS&tp,^,^) = - e -^(3i? 4 -4L> 2 ) = _^L( 3 £ 2 -^ 2 L 2 ) 
low 4w a 

= - a [3 - sin 2 (2$) (cos 2 >p sin 2 ^ 2 + sm 2 <p sin 2 ^ 3 ) - sin 4 i? sin 2 (2(^) sin 2 (^ 2 - v 3 )\ , (18) 

where L 2 = L 2 + L 2 + L 2 is the conserved squared angular momentum, and the first-order action unit a is given by 

a = enE 2 /4cj 5 . (19) 

The integral (|12fl . with the function given in (|18f) . has been integrated numerically and found to be identical with 
the corresponding five-dimensional integral over the 5-sphere S 5 , expressing AiS in terms of the five hyperspherical 
angles of the six-dimensional polar coordinates (cf. Appendix |B|) . However, even the four-dimensional integration 
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in (|12fl is more than needed. In fact, the same result is obtained if we replace the phase function under the CP 2 
integration (|12|) as follows: 

AiS(^) = -ff(3-cosV). (20) 
The integrals over z/2 and then become trivial and the result is found (with z — cos 2 cp) to be 

Mika/h) = 2 f 2 cos (p simp dip e- lka{3 - cos2lp)/R = [ dz e -*M3-*)/» = A [ e -ak*/K-i*/2 + e -m<r/h+in/2'} _ , 21 j 

Jo Jo fccrL J 

The replacement (I20|l is not obviously justifiable in a direct way, but a suitable reduction of the five-dimensional S 5 
integral yields exactly the integral i|21[l. as demonstrated in Appendix El Identically the same result is also obtained 
independently from the EBK quantisation of the perturbed Hamiltonian using Poisson summation and keeping only 
first-order terms in e; it is given in Eq. at the end of Sect. IIII CI 

The form (|21H of the modulation factor suggests a simple physical interpretation: the resulting two terms correspond 
to two separate families of orbits that survive the breaking of the U(3) symmetry. As we will see below, these are 
the circle and the diameter orbits with maximal and minimal (i.e., zero) angular momentum, respectively, at fixed 
energy. For k = 1 these are actually the shortest periodic orbits found in the perturbed system. In contrast to the 
four-fold degenerate unperturbed HO orbits, they only have a two-fold degeneracy, since only rotations about two 
of the Euler angles change their orientation. In general, the most degenerate periodic orbits in a three-dimensional 
system with spherical symmetry SO (3) can be rotated about three Euler angles. However, for each of the circle and 
diameter orbits one of the three rotations is redundant in the sense that it maps these orbits onto themselves: for the 
circles, it is the in-plane rotation around their centre, and for the diameters, it is the rotation about their own axis 
of motion. The loss of two degrees of degeneracy of these two orbit types with respect to the unperturbed HO orbits 
appears in the form of the factor h in their amplitudes in (|21|l , in agreement with the general counting of the powers 
of h in semiclassical amplitudes 0, . 

We will see in the next section that for values k > 3 and large enough values of e, there exist fully three-fold 
degenerate orbits, namely orbits with rational ratios to r : lo v = N : M of radial and angular frequencies with N:M> 
7: 3. They are created at bifurcations of the fc-the repetitions of the circle orbits with k = \M\ > 3. This happens, 
however, only at finite values of e, so that these orbits do not contribute in the limit e — > 0. 

Inserting (|21|l into ©, we find the following perturbed trace formula for the oscillating part of the level density: 



5g P ert(E) = 



k r 



k=l 

\k 



kS c 7r\ / kS d 37r 

cos — — + cos — — 

h 2 \ h 2 



c 



Aq, 2 f— lr 

" £ —r~ t sin (fcS c A) - sin(kS d /h) ] , (22) 



h 2 TT 



k=l 

where Sd and S c are the perturbed primitive actions of the diameters and the circle, respectively: 

S d = S - 3cr , S c = So - 2cr , S = 2ttE/lu . (23) 

These values will be confirmed from the general expressions derived in the next section. The perturbative trace 
formula (|22|l uniformly restores the oscillating part of the exact HO trace formula (JSJ in the limit e — * 0. In Section 
IIII El we will test its range of validity comparing against quantum- mechanical results. 

Note that - as is usual in perturbation theory - the semiclassical amplitudes and actions of the orbits contributing 
to (|22|l are correct only in the small-e limit. In general, one has to generalise the perturbative trace formula to a 
uniform version that, in the limit of large perturbation, goes over into the corresponding (extended) Gutzwiller trace 
formula. This uniformisation has been done for the breaking of U(l) symmetry in |l7j . and for some cases of U(2) 
symmetry breaking in [l^; one of the latter result applies, with suitable changes, also to some cases of SO(3) — > U(l) 
symmetry breaking. The symmetry breaking U(3) — > SO(3) under study here has not been treated in the literature 
so far. However, the simplicity of our above results makes the uniformisation particularly easy in the present case: 
since the modulation factor (|21|l already has its own asymptotic form - or, inversely speaking: since the asymptotic 
expansion of the integral in (|21|) in the limit of large a happens to be exact also for a — > - it will be sufficient to 
replace in (|22|l the perturbed actions Sd, S c and their semiclassical amplitudes by those valid for all values of e, which 
will be derived in the next section. 
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III. UNIFORM TRACE FORMULA FOR DIAMETER AND CIRCLE ORBITS 



In this section we will calculate the full actions and semiclassical amplitudes of the diameter and circle orbits of 
the Hamiltonian gjjg, valid for arbitrary values of e. The general trace formula for an arbitrary spherical potential 
in three dimensions has been given by Creagh and Littlejohn We choose a different approach here, which is 

in spirit that of Berry and Tabor , but goes beyond their leading-order approximation in taking into account the 
end-point corrections to a trace integral which yield precisely the circle and diameter orbit contributions. 

The classical Hamiltonian of the system (with mass m = 1) 

H{v,v) = E= l -p 2 + V{r) 1 V (r)= l -u 2 r 2 + l -er\ (24) 

is integrable due to the spherical symmetry of the potential, so that we can apply the standard EBK (or radial WKB) 
approximation to it. Writing H in terms of polar coordinates r = (r, 6, </>) and the associated canonical momenta 
p = (p r ,Pe,Pc/>), we have the usual form involving an effective potential V e //(r) that includes a centrifugal term: 

E = H{r lPri L) = \ P 2 r + V eff (r) , V eff (r) = V(r) + ^ , (25) 

where p r is the radial momentum. The three independent (and Poisson-commuting) constants of the motion are the 
energy E, the total angular momentum L 2 , and its z component L z = p^. The momentum pg can be expressed in 
terms of the latter two as 

pe = \j 'L 2 - L 2 Jsm 2 9 . (26) 

Before we specialise to our particular potential V(r) in (|24|l . we briefly recall the radial EBK quantisation and derive 
from it a general trace formula for an arbitrary spherical potential, starting from the Hamiltonian (|25|l . 



A. EBK quantisation of a spherical system 

In the standard radial EBK method 0, 13 > one quantises the three following action integrals 

S r = <j>p r dr = 2nh(n r + 1/2) , n r = 0,1,2,... (27) 

h = ^ <£p g d6 = h{n + 1/2), n g = 0,1,2,... (28) 

J = L z = hm, m = 0,±1,±2, . . . (29) 

where the Maslov index 1/2 in (|27J) is correct only for smooth potentials. Since pg in l|26(l does not depend on the 
potential, the integral for Ig can be done once for all and yields 

Ig = L-\L z \>0, (30) 

so that the quantisation condition for L is given by 

L = H(n e + \m\ + 1/2) =H(£ + l/2), I = n e + \m\ = 0, 1, 2, . . . (31) 

in terms of a single angular momentum quantum number I. The relation Q31JI includes the so-called Langer correction; 
the quantised squared angular momentum L 2 = h 2 (I + 1/2) 2 = h 2 (£ 2 + £ + 1/4) agrees with the exact quantum- 
mechanical value h 2 1(1 +1) in the limit of large i. Solving 1251) for p r , we can write the radial action as 

S r (E, L) = <j) dr^2E ~2V{r) - L 2 /r 2 = 2ith{n r + 1/2) , (32) 

showing that the quantised energies will only depend on the radial and angular momentum quantum numbers n r and 
£; they have, of course, the usual m-degeneracy di — (2£ + 1), since —£< m < +£. 
Inverting the relation (|32J) . we may rewrite the Hamiltonian Ij25(l in the form 

E=H(S r ,L). (33) 

Inserting the right-hand side of H27J1 and (|31Jl into (|33|l . we obtain the EBK-quantised eigenenergies: 

E^ b r k e =H(2Trh{n r + l/2),h(£+l/2)), I, n r = 0, 1, 2, . . . (34) 

They can in general only be obtained by numerical iteration after doing the radial action integral l|32|) over r within 
the classical turning points. 
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B. Introduction of scaled variables 



Before continuing, we simplify the situation by a scaling of the energy. In principle we have to vary the three 
parameters E, uj, and e to study the dynamics of our present system. However, we can introduce a scaling of the 
energy in such a way that the classical dynamics only depends on one single parameter, a dimensionless scaled energy 
e. If we multiply by the factor e/w 4 , we can write the r.h.s. in terms of scaled coordinates and momenta pi and 
a scaled time r: 

"=^-Xi, Pi = —rPi, r = wt, (35) 



so that the scaled energy e becomes 



c 



1 _o , , 1 .o , , I 2 , . 1 o 1 



e = ^E=-p 2 + v(q) = -q 2 + v(q) + -^, v(q) = - q 2 + - q 4 . (36) 
where the dimensionless scaled angular momentum I is given by 

UJ 3 

I -LI a, s= — . (37) 
e 

In 13tj|) q is the scaled radial variable and q — p q the scaled radial momentum, and the dot means the derivative with 
respect to the scaled time r. s is the action unit. 

This brings about a considerable simplification of the classical dynamics: we only need to vary one parameter e; at 
the end of our calculations we just have to remember that energies are measured in units of w 4 /e, angular momentum 
and actions in units of s = w 3 /e, times in units of l/u>, etc. In the following we shall give all quantities as functions 
of the dimensionless scaled variables e and I. (Other scaled dimensionless quantities such as actions, periods etc. will 
be denoted by lower-case letters.) 

The scaling in l|3()|) is specific for our present Hamiltonian l|24f) . It can, however, easily be modified for HO potentials 
perturbed by arbitrary central potentials which are pure power laws in q. The results of this subsection can therefore 
easily be generalized to the corresponding suitably scaled potentials v(q). 

C. Density of states for an arbitrary shperical potential 

We now take the spectrum (|34H as a starting point to write down the density of states in the EBK approximation: 

oo oc 

9ebk(e) = Y, £( 2 * + !) % - <5) • (38) 

n r =0 e=o 

Next we apply Poisson summation |25j | to convert the sums over n r and I into integrals: 

9ebk(e)= E / ^(2^+1)/ dn r 6(e-e e n b r k e )e i2 ^ N ^ +M ^ + ... (39) 

AT 1, JO JO 

N — — oc M— — oc 

Due to the finite lower limits of the summations in (|38|l . there are boundary corrections to (|39|l . which we have 
indicated by the dots. We shall comment on them in a moment. Using (|27|l and (|3ip. we substitute the variables £ 
and n r by the classical actions L and S r , using d£ — dL/h and dn r = dS r /2Trh. Then g e &fc(e) becomes, expressed in 
terms of the dimensionless scaled variables, 



ebfc (e) = -i-^V V / dl2l d S ^(e-^( Sr ,0)e-^' +2 ^]/ R -^ +M ) + ... (40) 



Hereby s r — S r /s is the scaled radial action integral and h (s r , I) the scaled Hamiltonian (|33|l . Note that the integration 
limits have been imposed by the energy conservation; / m (e) is the maximum scaled angular momentum at fixed energy. 
The lower limits are, strictly speaking, hs/2 for the integral over / and tms for the integral over s r . Replacing them 
by zero corresponds to neglecting corrections of higher order in h. In the following, we keep terms up to order H with 
respect to the leading factor oc fir 3 in (|40|l . There are exactly two corrections of relative order h in what is indicated 
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by the dots above: one is a boundary correction from n r = to (|39|l . and the second is coming from the lower limit 
■whs of the s r integral in l|4UI) . A short calculation shows that these two corrections cancel identically. All other terms 
neglected in H4()[) correspond to corrections of relative order ft 2 or higher. 
We now use the relations 

S(e-h(s r , I)) = \s' r (e, l)\S{s r - s r (e, I)) (41) 

and 

s' r (e,l) = ^^-=t r (e,l). (42) 

Here f r (e,/) = u>T r (e,l) is the scaled period of the classical motion at fixed values of e and I, which is the energy 
derivative of the corresponding scaled action integral s r (e, I): 

'1-2 



p q {e,l)dq = 2 dq^2e - 2 



s r (e,l) = d>p q (e,l)dq = 2 / dq^/2e - 2v(q) - l 2 /q 2 , (43) 



tr (e,l) = = 2 rdq - 1 (44) 

K de J qi J2e - 2v(q) - l 2 /q 2 

Here q\ and q2 are the scaled lower and upper turning points of the radial motion, respectively, which both depend on 
e and I. The above integrals can in many cases be expressed in terms of complete elliptic integrals. (Exceptions are 
the HO and Coulomb potentials where they become simple algebraic functions of e and /.) For our present potential 
v(q) in (J2EI, the analytical expressions for ll4^|l . ll4^|l and other quantities of interest are given in the Appendix El 
One important result derived there is the Taylor expansion of s r (e,l) around 1 = 0, whose first terms are 

s r (e, = s r (e, 0) - 7T I + a(e) I 2 + 0(l 3 ) , (45) 

where a(e) is given in i|C9(l . The structure of (|45|l - but not the explicit form of a(e) - appears to be a general result 
valid for all regular central potentials v(q) with a minimum at q = (and hence not for the Coulomb potential, see 
Appendix ID 2|l . but we were not able to prove it in the general case. The relevance of the linear term — ttI in (|45() will 
become clear in the next subsection. 

Using the above relations we can now do the integral over s r in 140J) exactly, due to the delta function, and obtain 
the following "EBK trace integral" : 



1 5 00 00 pi (e) 

l!ihk { e ) = J- < £- £ J2 (-l) N+M dllt r (e,l)e^ NS ^+^m + o{h^). (46) 

n £ AT=-oo Ai=-oo ^° 



The phase of the integrand of (|4l)|) can be interpreted as the full action Snm{&, = NS r (e, I) + 2nAIsl, divided by 
h, of a given classical orbit labelled by M and N, consisting of the radial part NS r (e, I) and the angular part 2nMsl. 
When doing the full double summation over all M and iV and the integration over / exactly, (I4(j|l yields the EBK 
spectrum (|34() to leading order in h. Note that the limits of the / integral are the cases of zero angular momentum, 
which corresponds to the diameter orbits, and its maximum value l m (e) at a given energy, which corresponds to the 
circle orbits. The latter have the radius go at which the effective scaled potential v e ff(q) has its minimum; for this 
motion the radial action is zero: S r (e,l m (e)) = 0, since all energy is in the angular motion. All contributions with 

< I < l rn (e) to the integral correspond to motion which has both radial and angular components. 

The formula (|46|l . valid for arbitrary (but correctly scaled) spherical potentials v(q) with smooth walls, is in principle 
a trace formula, but it does not yet have its characteristic form. That form is obtained by evaluating the integral over 

1 in the semiclassical limit h — > by evaluating its leading contributions from the critical points of the phase function 
Snm (e, in the exponent of the integrand. To leading order in ft, these are stationary points - as far as they exist. 
Using the standard stationary-phase evaluation of the integral around the stationary points, one obtains contributions 
to g e bk(e) with amplitudes proportional to ft -5 / 2 , as expected for the fully three-fold degenerate orbits in a spherical 
system (cf. the discussion at the end of the previous section). Next to leading order in ft, one obtains contributions 
from the end points of the integral (see, e.g., j26j), sometimes referred to as "edge corrections". As already announced 
above and shown explicitly below, these correspond here to the diameter and circular orbits, giving contributions with 
amplitudes proportional to ft~ 2 . 

A special contribution to g e bk(e) comes from M = N = 0. As generally proved by Berry and Mount |2jJ, this must 
be the Thomas-Fermi (TF) value of the density of states, which is the leading contribution to its smooth part. From 
(gBJ we obtain with M = N = 

9%l(e) = ^^f m{e) dl 2 t r (e,l). (47) 
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The integral of (|44|) over I 2 is straightforward. Note that the contribution from the upper limit l m (e) in (|47|l is zero, 
since the turning points coincide: <?i(e, l m ) = 92(6, l m ) — Qo- The result is 



<l 



2uf 
irh 3 e 



2 dqy/2e~2v{q)=g TF (E). 



(48) 



Here q§ is the upper turning point for I — 0. (The lower one is zero: q® = 0.) That l|48|) really is equal to the TF 
density of states follows from its general definition 



grF(E) = 



1 



d A p I (Tr S(E - ff(r,p)) = 



(2tt/i) 3 



d 3 p A*(e-/>(q.p)). 



(49) 



Using polar coordinates for q and p and doing the p integration leads to H48f) . The analytical expression valid for 
the potential is given in (|C16|) . In the introduction we have stated that the average part of the density of states 
is generally given by the ETF approximation which contains h corrections to the TF limit. For the spherical three- 
dimensional HO, these corrections are of relative order h 2 , as seen in (0. For the present perturbed HO potential @ 
we will see that the TF approximation is sufficient to reproduce the average part of the quantum-mechanical density 
of states, at least up to the energies for which the quantum spectrum is numerically available. 
The stationary condition for the phase function in (j46(l at a point I = Inm reads 



^-[Ns r {e,l) + 2TrMl] 
01 



= N 



ds r (e, I) 



01 



2ttM = . 



Due to energy conservation we may write 



~ Oh 
de — dh (s r , I) = — — 
0s r 



dh 

ds r + — 
al 



dl = 0, 



so that 



dsr(e,l) to<t>(e,l) 

= — Z7T ■ 



dl u r (e, I) ' 

where the frequencies of angular and radial motion arc defined as usual by 



dH(S r ,L) 
dh 



Lo r (e, I) = 2ir 



dH{S r ,L) 
dS r 



With this, the stationary condition H5()[) becomes 

w<p(e, Inm) 



T r (e,l NM ) _ M_ 
uj r (e,l NM ) T^e, Inm) N 



(50) 



(51) 



(52) 



(53) 



(54) 



which is the periodicity condition for the "rational tori" as the most degenerate classical orbits || • Clearly, N and M 
must have the same sign, since frequencies and periods are positive quantities. Because of (I45|l . the lower integration 
limit I = in (|46|l is a stationary point corresponding to the diameter orbit with N : M = 2:1. However, the 
stationary-phase integration of (|46|l for I = I21 = yields a zero contribution due to the factor I in the integrand, so 
that the only semiclassical contribution due to the diameter orbits is the end-point correction for I — discussed in 
the next section. As we shall see further below, solutions of (|54|l with < Inm < lm(e) do not exist for all energies e 
and for all pairs N, M of integers. Therefore we postpone the contributions of the rational tori to (|46|l and concentrate 
first on the diameter and circle orbits. 

Before developing the corresponding trace formula for our present system, valid to all orders in e, we want to 
establish here the connection to the first-order perturbative approach used in Sect. [H] From the Taylor expansions 
given in Appendix El we obtain for the radial action integral the following first terms: 



2SJE,L) =2n{ L 

UJ 



4w 5 



(3E 2 



2 T 2\ 



lu z L 



0(e 2 



(55) 



In the second term we recognise precisely the first-order action shift given in (|18|l . Inserting into it L = and the 
leading-order expression for L m — E jio given in l|C22|) . we obtain the first-order action shifts of the diameter and circle 
orbits, respectively, given in l|23|l . We now insert l|55|l into the expression l|46|l for the density of states, neglecting the 



10 



terms of higher order in e and keeping only the zero-order terms in the amplitude of T r — tt/ui, given in (|C18() . and 
of the upper integration limit L m — E/to. Noting furthermore that for the unperturbed HO orbits the ratio T r /T^, 
becomes equal to 2, so that the resonance condition 1)54(1 implies N = 2M, we assume for the moment that all other 
combinations of N and M in this limit may be neglected. Writing N = 2M = 2k, we obtain the following result for 
the density of states up to first order in e (note that the linear terms in L cancel in the exponent of the integrand) : 

ff W( jB) = _L_ V (-l) fe e i2fffeE /^ / dL 2 e -ik*(3-u 2 LVE*w^ (56) 

where a is the quantity defined in (|19fl . Using the substitution z — (loL/E) 2 , the above integral becomes identically 
equal to the integral for the perturbative modulation factor ((21(1 . and the result 1(56(1 is precisely the perturbed density 
of states defined in with the oscillating part given in l(22ll . Rather than proving at this stage that all contributions 
to 146(1 with N 2M either cancel or are negligible to the leading order in h, we now go on to derive the full trace 
formula for the diameter and circle orbits, doing the summations over all N and M exactly in the semiclassical limit. 



D. Semiclassical trace formula for diameter and circle orbits in the perturbed harmonic-oscillator potential 

Equipped with the above results, we are now in a position to evaluate the full trace formula for the contribution 
of diameter and circle orbits for our present Hamiltonian @. The smooth TF part of the level density, qtf(E), is 
given explicitly in the appendix [U] The contributions of the rational tori, which bifurcate from the circle orbits only 
for high enough energy and repetition numbers, will be discussed in Section IIII Fl We therefore now evaluate the 
end-point contributions to the integral in l(46|) in the semiclassical limit h — ► 0. 



1. Diameter orbit 
The lower end point I — in 14611 yields the asymptotic contribution 



5gd(e) 



2neh 2 



T,.(e,0) £ (-l) N e 



i[NS r (e,0)/h+TT/2] 



(57) 



N=-oo 

where we have left out the contribution from N = which contributes to the TF part, and defined the quantity 

(-1) M 



uiv(e) = lim 



' X] I A - Is 



.M= — oo 



(58) 



Using the identity |28 



E 



(-i) 



M 



we find 



sin(z) ^— ' (z — Mit) 



Uw(e) = lim 
i^o 



(z ^ tvk, n£Z) 



(N ds r \ 



We now exploit the result l(45|) from which we find 



ds 



^ = -TT + 2a{e)l + 0{l 2 ). 
al 



(59) 



(60) 



(61) 



For odd values of N, the sin function in ((60(1 gives always a nonzero denominator in the limit I = 0, so that itAr(e) 
becomes zero. For even N = 2k we get in the limit I — > 



mii I ) = sin(-fc7r + 2kla(e)) = (-l) fc sin(2fcZo(e)) — ► (-l) fc 2Jfc/»(« ) 



(62) 
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so that we obtain 

U ^) = ^y («) 

Inserting this into (|57|l using (|C21|) . leaving out the contribution k = which is contained in <?TF(e), we obtain the 
final contribution of the diameter orbits to the density of states: 

^(e) = - i ^E^ L -(^(e)A), (64) 
v ; fe=i 

where the actions Sd(e) = 2S r (e, 0) are given explicitly in (|C12(1 . In the low-energy limit we obtain with (|C9|) exactly 
the diameter contribution to the perturbative trace formula l|22|) . 

Note that the presence of the linear term —irl in the quantity 1|45[) is instrumental in annihilating the contributions 
with odd N to the sum in (|57|l and hence in establishing the fact that only an even number of radial oscillations yields 
a physical periodic orbit with angular momentum L = 0. 



2. Circle orbit 



The upper end point I — Z m (e) in (|46[) yields the asymptotic contribution 



9c (e) = ±L m (e)T r (e,l m ) £ (-i^P^e)/^] £ -Lii ^. (65) 



Using ISHI and noting that T^(e,l m ) = T c (e), we get 



ds r 
~dl 



2iTT r (e,l m ) 

'^UeT- (66) 



Employing l|59|) again to perform the N summation, writing \M\ = k and omitting the k = term, we obtain after 
some manipulations the contribution of the circle orbits to the density of states 

= T c (e)L . (e) g (-1)* (e)/ft) (67) 

where S c (e) is given explicitly in l|C22l) . and the quantity w(e) is defined by 

w{e) = = — f—r^;- (68) 

±r{e,lm) w</,(e,t m j 



In the limit e — > this quantity is found with the r.h.s. of the analytical results l|C23|) and (|C25|I to go like 

w(e) = 2(1 + e/4+...) =2+ E+... (69) 

Expanding the denominator in (|67|) . using the above w(e), up to first order in e brings the amplitude (|67|l exactly 
into that of the circle orbit contribution to the perturbative result <|22|) • 

Note that the denominator under the sum in (|67|l looks exactly like that in the Gutzwiller trace formula 0, |2j| 
for an isolated stable orbit, whereby irw(e) corresponds to one-half of the stability angle. Indeed, the semiclassical 
amplitude in (|67|l diverges when kw{e) becomes an integer n£Z. At the corresponding nonzero energies, the circle 
orbit bifurcates; the condition for this to happen is exactly the resonance condition (|54|l for the rational tori with 
M = k, taken at Inm = ^m(e). (We come back to this point in Section Ull FJ ) Using 1C23|) and (|C25fl . we find that 
w(e) is restricted to the following range: 

2 < w(e) = Tc ^ < V6 for < e < oo . (70) 

The smallest k for which we can have 2 < N/k < ^6 = 2.4494897 with N <E N is k = 3 with N = 7, so that 
w(er:3) = 7/3, which happens at the scaled energy e7 : 3 ~ 7.670 (see Table |l| below). This means that the shortest 
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new orbit is the 7:3 torus, bifurcating from the 3rd repetition of the circle orbit. The next bifurcation (from k = 4) is 
that of the 9:4 torus at eg : 4 ~ 2.0967. Still longer orbits bifurcate at lower energies, but from higher repetitions of the 
circle orbit. Therefore, we can ignore the bifurcations at low energies and small repetition numbers k and concentrate 
on the contributions of the diameter and circle orbits to the density of states. 

Let us finally write down the trace formula which we have obtained for the oscillating part of the density of states: 



5g(e) ~ H*(«0 sm(kS c (e)/H) + A d k (e) sm(kS d (e)/H)] , (71) 
fc=i 

where the semiclassical amplitudes are given by 

(e) = T c (e)L m (e) (-1)" T d (e)^ (-1)^ 

k[ ' ttH 2 sm[kirw(e)} ' k[ ' Aneh?a{e) k ' 1 ' 

In the following we shall test this formula against the quantum-mechanical density of states. As we have shown above, 
(I71|l goes over into the perturbative trace formula llL'l'l) in the limit e — > 0. This confirms the assumption, made at the 
end of Sect. IIIlUl that the contributions with N ^ 2M to (|46|l do not contribute in this limit. 



E. Numerical tests of the trace formulae 



In this section we compare our semiclassical results with those obtained from the exact quantum-mechanical spec- 
trum (obtained numerically by solving the radial Schrddinger equation on a discrete mesh). In order to focus on the 
coarse-grained shell structure, we convolute both results over the energy E with a normalised Gaussian of width 7. 
The quantum-mechanical density of states Q then becomes 

9^E) = ^J=^- {E - En)2h2 - ( 73 ) 

7V7T ' 
v n 

In order to obtain its oscillating part, we subtract from it the TF expression gTF{E) which we can calculate analytically 
(see Appendix EJl ■ The semiclassical trace formula (J2J becomes, after doing the convolution in stationary-phase 
approximation (cf. 0, Sec. 5.5) 




-200 I , , , , J 

175 180 185 ^ „ n 190 195 200 

E [huj] 

FIG. 1: Density of states for the Hamiltonian (J^J with e = 0.01, Gaussian-averaged with a width 7 = 0.5hu, versus 
energy E (units: hui). Dashed lines (red): quantum-mechanical results. Solid lines (black): semiclassical results using 
the uniform trace formula 17H . 
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Sg sc (E) ~ J2Ao(E) e -hT po (E)/2H? CQS [Spo{E)/h _ CTpo7r/2] 5 (74) 

po 

so that orbits with longer periods T po will be exponentially suppressed. 

In Fig. ^ we show a comparison of the results obtained with the uniform trace formula (|71|l for the case e = 0.01 
with those obtained from the exact numerical quantum spectrum. The width of the Gaussian smoothing was chosen 
to be 7 = 0.5Hu>; in the semiclassical result the harmonics k > 2 then do not contribute noticeably. The agreement is 
seen to be perfect; tiny differences can only be noted near the beat minima on the amplified scale below. 

In Fig. |21 we show a similar comparison for 7 = O.lhui, exhibiting much fine structures. We now start to resolve 
the spectrum with a higher resolution; note that we still only plot the oscillating part Sg(E) of the density of states. 
In the top panel, the semiclassical result includes harmonics (i.e., repetitions of the two orbits) with k < k m = 7. 
This is evidently sufficient to reproduce the quantum-mechanical oscillations up to E ~ Abhw with a high accuracy. 
In the centre panel, we have added two more harmonics, going up to k < k m — 9. Here we recognise the appearance 
of divergences at the scaled energies e = 0.361712 and e = 0.437867 which correspond to the bifurcations of the 19:9 
and 17:8 tori from the repetitions of the circle orbit with k = 9 and k = 8, respectively (cf. Table [I] below) . The 
small wiggles close to the divergences, with decreasing amplitudes when going away from them, are due to the missing 
contributions from those torus orbits in the trace formula (|71|l (see the following section). 




10 20 ^ r , n 30 40 50 

E [ha;] 




33 34 35 36 37 38 39 40 41 42 43 44 45 46 

EM 

FIG. 2: Same as Fig. Q but with 7 = O.lhuj. Top panel: including up to k m ~ 7 harmonics in the trace formula l|71|l . 
Centre panel: adding harmonics up to k m — 9. Note the divergences due to the bifurcations of the 17:8 and 19:9 to ri a t 
the energies E ~ 43.79 and E ~ 36.17, respectively. Bottom panel: result obtained from the EBK trace integral l|46|l . 
doing the integration over I numerically, with — k m < M < +k m ., —2k m < N < +2k m using k m = 9. 
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In order to anticipate the results of a suitable uniform trace formula including the contributions of the tori, we show 
in the bottom panel of Fig. [2] the result of the EBK trace integral (|46() obtained by a numerical integration over the 
variable I. Hereby the summations over N and M have been done in the limits — k m < M < +k m , —2k m < N < +2k m 
using k m = 9. [The period to be used in the Gaussian damping factor of (|74f) is in this case given by T po — NT r (e, L).] 
Now the divergences and all other discrepancies have disappeared; the difference to the quantum-mechanical result 
can hardly be recognised. The integral l|46l) therefore is also a uniform expression for the semiclassical density of 
states, containing the contributions of all periodic orbits. The numerical evaluation becomes, however, very time 
consuming for even finer energy resolutions (i.e., still smaller values of 7). 

Our main result here is that up to a rather fine resolution, the shell structure of the present perturbed HO 
system is dominated by the two-fold degenerate families of diameter and circle orbits, although these are not the 
most degenerate orbits. This is a rather unusual situation, due to the fact that the shortest tori with the highest 
degeneracy are considerably longer (by a factor > 10) than the primitive diameter and circle orbits. 

Let us finally test the range of validity of the perturbative trace formula ll'L'l) . In Fig. [3] we compare its results to 
that of the full trace formula l|71|l for the two values e = 0.001 (upper panel) and e = 0.01 (lower panel). As expected, 
the two curves agree in the limit e — > 0. However, already for e = 0.001 the amplitude modulation of the perturbative 
result deviates so much that the position of the first beat node is predicted by about 10% too low in energy. For 
e = 0.01, the perturbative result becomes so bad that it predicts the second beat node approximately at the position 
of the correct first one. These results underline the fact that in [2l|, the quantum results for weak interactions could 
be qualitatively reproduced by the perturbative trace formula, but not quantitatively with correct positions of the 
beat nodes. A more detailed analysis of the selfconsistent HF results of [21| using the uniform trace formula l|71|) is 
in progress |30j . 



400 





FIG. 3: Same quantity as shown in Fig.Q Here we compare the results of the perturbative trace formula 1221 (dashed 
lines, red) and those of the uniform trace formula 17H (full lines, black), Gaussian-averaged with 7 = 0.5huj. Upper 
panel: for e = 0.001; lower panel: for e = 0.01. 



F. Bifurcations of rational tori from the circle orbit 



As we have seen in the previous section, periodic orbit with both angular and radial motion appear as the rational 
tori fulfilling the periodicity condition 154|) . For the Hamiltonian this condition can only be fulfilled in the limits 
given in JTUJ). This means that stationary points Lnm of the L integral in l|46[l satisfying i|50|) and hence l|54(l only 
exist if 2 < \N : M\ < ^/6. Hereby \M\ corresponds to the repetition number k of the circle orbit that undergoes a 
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bifurcation exactly when Lnm is equal to the maximum angular momentum L m . In Tabled we list some of the lowest 
energies e^-M at which this happens for the smallest integers M = k. Each of the N:M tori only exists above its 
bifurcation energy, i.e., at energies e > cn-.m- 



M = k 


N 


e-N-.M 


Inm 


,(e=oo) 
NM 


9 


19 


0.3617121 


0.3370 


0.276 


8 


17 


0.4378671 


0.4032 


0.316 


7 


15 


0.5527726 


0.5009 


0.360 


6 


13 


0.7441150 


0.6584 


0.420 


5 


11 


1.1174198 


0.9508 


0.504 


4 


9 


2.0966668 


1.6553 


0.633 


3 


7 


7.6699999 


4.9331 


0.864 



TABLE I: Bifurcation energies e^-.M and angular 
momenta Zjvj\/ of the ratio nal N :M t ori at which 
the stationary conditions 1501 and 1541 1 are ful- 
filled for the lowest values of M and N. In the 
last column we give the angular momenta Inm for 
the limiting case e — > oo in which t hey do not 
depend on the energy (cf. Sect. llllTjl . 



The contributions of the bifurcated tori to the semiclassical trace formula is obtained by a stationary-phase evalu- 
ation of the I integral in Ij46(l , leading to 



(e)= £ ^ M ( e )(-l) w+M cos(%M + j). (75) 

AT>2A/>0 V " 



The summation goes only over those pairs N, M of positive integers for which the stationary condition l|54|) can be 
fulfilled. (Note that the contribution of the corresponding pairs of negative integers is included in the amplitudes by 
a factor 2.) The actions of the tori are given by 



Snm (e) = NS r (e, l NM ) + M 2nL NM (e) , (76) 



and their amplitudes by 



2\/2s I 



d 2 s r (e, I) 



dl 2 




(77) 



where s is the action unit given in l|37|l . The Heaviside step function in l|77|l is defined as 

6{x) = for x < , 6{x) = 1/2 for x = , 6{x) = 1 for x > , (78) 

where the value 1/2 for x = accounts for the fact that one obtains only half a Fresnel integral when the stationary 
point is the upper end point of the integral in 146|) . The power hr^l 2 in the amplitude Anm{g) is characteristic of 
the three-fold degenerate orbits in a three-dimensional spherical system 0, 0] ■ 

Eq. Ij75(l corresponds to the standard Berry- Tabor trace formula for the most degenerate orbits in an integrable 
system. Here, however, each contribution to (|75|l is only valid if the energy e is sufficiently larger than the corre- 
sponding bifurcation energy ejv:M- In the neighbourhood of the bifurcation energies, we have to replace the sum of 
the diverging circle orbit contribution and the corresponding torus contribution with \M\ — r by a common uniform 
contribution. As it is well-known from the semiclassical theory of bifurcations (l6Hl!^| . one has to include hereby also 
the contributions of so-called "ghost orbits" , which are the imaginary (or complex) continuations of the bifurcated 
orbits (here: the tori) to the other side of the bifurcation (here: e < cn-.m)- 

It is, however, not the purpose of our present paper to discuss in more detail the uniform approximation suitable 
for this situation. The present scenario is actually almost identical to that described in [3ll | for the bifurcations of 
tori in the integrable Henon-Heiles system. Although there the orbit from which the tori bifurcate is an isolated one, 
while it here is a family of circle orbits, the uniform trace formula given in |3l| can be applied to the present system 
in a straightforward manner. We refer the interested reader to this paper for all the detailed formulae as well as for a 
numerical demonstration showing how diverging discrepancies of the type seen in the centre part of Fig. [21 disappear 
when a proper uniform trace formula is used. Its result would here be practically indistinguishable from the result 
shown in the bottom part of Fig. |2 in which we have evaluated the integral in (|46|l numerically instead of using the 
stationary phase approximation with end-point corrections. 
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G. The limit of the purely quartic oscillator 

In this section we discuss briefly a purely quartic oscillator potential, i.e., we start from the Hamiltonian 

H(r,p) = l -p 2 + ( -r 4 = E. (79) 

Here both the energy and e can be scaled away. After dividing the above equation by E and introducing the scaled 
variables 

CH = (|) 1/4 r, , r = (eE)V* t, 1 = e^E^L , (80) 
we obtain the "unit energy" equation 

l = \ e + \ qi + ^' (81) 

where the dot again means derivative with respect to the scaled time r. Hence the classical equations of motion 
become independent of energy and e. After obtaining all the classical results, we can reintroduce E and e just by 
remembering that lengths scale with (E/e) 1 ^ 4 , momenta with y/E, times (periods) with (_Ee) -1 / 4 , actions and angular 
momenta with e~ 1 / 4 E 3 / 4 . 

We obtain the semiclassical trace formula from the results of the previous sections and of the appendix O simply by 
taking the limit e — > oo, which means practically by extracting the leading terms for e — > oo. In this way, we obtain 
the following primitive periods and actions for the circle orbit: 

T C (E) = 2 7 re- 1 / 4 (4£/3r 1 / 4 , S C (E) = 2nL m (E) = 2 7 r e - 1 / 4 (4£/3) 3 / 4 , (82) 

and for the diameter orbit: 

T d {E) = 2V2e- 1/4 K Q E- 1/4 , S d (E) = 2y/2 e-^ 4 (A/3)K Q E 3 / 4 . (83) 

Here Kq is the constant elliptic integral 

K = K(k ) = 1.8540746773. . . , k = 1/V2. (84) 

Furthermore we get 

T r (E,L m )=ir(3eE)- 1 / 4 . (85) 
The resonance condition (|54|l becomes independent of the energy 

T<p(E,l NM ) n > N 

UE^j =W ^=M> (86) 

since both periods scale with the same power of the energy. For the frequency ratio we find the limits 

2 < w(l NM ) = N:M < V6 for < l NM < l m = (4/3) 3/4 = 1.240806 .. . (87) 

Note that the upper limit l m never becomes a stationary point, so that no bifurcations of the circle orbits occur in 
this system. The values of the stationary points < Inm < lm (in scaled dimensionless units) for the lowest values of 
N and M are given in the rightmost column of Table [U as we just have seen, they do not depend on the energy E. 

The result w(l m ) — y/6 is a special case of the general result, given in [32j . that for a potential V(r) = Uq r@ one 
has w(l m ) = y/P+2. 

For the amplitudes of the circle and diameter orbits in the trace formula (|71|l . which keeps the same form, we obtain 
here 



h 2 V3 V e sin(fc7rV6) ' ^ 2 tt 2 V e k 

The actions S c and Sd to be used in l|71|l are those given in (|82|l and 1)83(1 . Finally, the TF level density is found to be 

zK/2 E 5 / 4 
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The tori, which in the perturbed HO system discussed in the last section occur only at sufficiently high energies 
after their bifurcations from the repeated circle orbits, exist in the present system at all energies. However, they affect 
again only the finer quantum structures of the density of states since they only exist for N : M > 7:3 due to the 
selection by their resonance condition (|87|) (cf. also Table HJ). This is demonstrated in the following two figures. 

In Fig. we compare the results of the semiclassical trace formula l|71|) for the diameters and circles using the 
above quantities to the exact quantum-mechanical density of states, both Gaussian-averaged over a width 7 = 0.5 (in 
energy units). The agreement is again perfect even on the enlarged scale in the lower part of the figure. The reason is 
that for the chosen energy resolution with 7 = 0.5, only repetition numbers \M\ — k < 3 contribute, for which there 
exist no tori (see Table HJ. 
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FIG. 4: Density of states for the Hamiltonian l|79|l with e = 0.01, Gaussian-averaged with a width 7 = 0.5 energy units. 
Dashed lines (red): quantum-mechanical results. Solid lines (black): semiclassical results using the uniform trace formula 
1711 and the actions and amplitudes given in 1821 . 1831 and 1881 with k < 2. 
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FIG. 5: Same as Fig. [3] but with 7 = 0.2 energy units. Repetitions of the diameter and circle orbits with k < 5 are 
included in the semiclassical result. 



Fig. [3] shows the same comparison for 7 = 0.2. Here some small discrepancies appear which become more visible 
with increasing energy. They are due to the missing contributions from the tori with \M\ — k > 3. Since there are no 
bifurcations involved in the present system, the contributions of the tori are given by the Berry- Tabor trace formula 
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(|75|l . employing the appropriate expressions for the actions and periods used therein. The analytical expression of 
the action integral S r (e,l) in terms of elliptic integrals is the same as for the perturbed HO potential, as given in 

Eqs. (|Cip - l)C7j) . except that here the term q 2 /2 in the potential and the corresponding contribution s"r\e,l) must 
be omitted. 

Even at the rather fine resolution obtained in Fig. \5\ our main result is the observation that the circle and diameter 
orbits dominate the shell structure also for the purely quartic oscillator. 

IV. SUMMARY AND CONCLUSIONS 

Motivated by recent Hartree-Fock (HF) calculations of harmonically trapped fermionic gases with a repulsive 
interaction |2l| , in which pronounced super-shell effects were obtained, we have developed a semiclassical trace formula 
for the density of states of the anharmonically perturbed three-dimensional isotropic harmonic oscillator (HO) defined 
in We find that its gross-shell structure is dominated by the periodic orbits of diameter and circle shape, whose 
interference leads to super-shell beats which explain the findings from the numerical HF results. 

In a first step, we have used the perturbative approach of Creagh to describe the symmetry breaking U(3) — > 
SO (3) for weak anharmonicities e, resulting in the perturbative trace formula 1|22|) . It uniformly restores the U(3) 
limit, yielding the exact trace formula (jSJ) of the unperturbed HO system in the limit e — > 0. For the derivation of the 
perturbative result l)22[l. one must in principle integrate over the 4- fold degenerate families of unperturbed HO orbits 
which cover the manifold CP 2 . As shown in App.[51 however, it turns out to be easier to integrate over the energy 
shell in phase space, which is a five-sphere S 5 , whereby the integration can be reduced, under exploitation of the SO (3) 
symmetry, to the one-dimensional integral given analytically in 1|21|) . An interesting result hereby is the fact that its 
two contributions, corresponding to the diameter and circle orbits, already have the characteristic form occurring in 
the general trace formula @. Usually, this form is obtained from a perturbative trace formula only asymptotically 
in the semiclassical limit h — > 0; this happens, e.g., also when one investigates the corresponding Hamiltonian @ in 
two dimensions or other perturbed 2-dimensional systems [|| 12^, . 

Next we have used the standard EBK quantisation and the Poisson summation formula to derive a general semiclas- 
sical trace formula for arbitrary central potentials in terms of a one-dimensional integral over the angular momentum 
L of the classical orbits, given in Q46[l. Exact integration and summation over all M,N yields, to leading order in h, 
the EBK spectrum. Its asymptotic evaluation in the limit H — > yields the standard-type contributions to the trace 
formula The end-point contributions yield the diameter orbits with L = and the circle orbits with maximum 
angular momentum; both contributions are of next-to-leading order in h, whereas the leading-order contributions 
corresponding to the typical rational tori of integrable systems according to the Berry- Tabor theory [8j come from 
the stationary points inside the integration interval. One interesting mathematical aspect is the mechanism by which 
the ratio cu r : lo^ = 2 of radial to angular frequency of the diameter orbit is naturally selected by a property of the 
general radial action integral S r (E, L) — § p r (E, L, r) dr as a function of angular momentum L which appears to be 
universal for regular spherical potentials. It holds also for a particle in a spherical box with specular reflection, for 
which we have re-derived in Appendix ID II the trace formula given by Balian and Bloch Q from our general formula 
(|46H . taking into account the changes of the Maslov indices due to hard- wall reflections). The Coulomb potential, 
for which u> r : uj^ = 1, is discussed in Appendix ID 21 where we also re-derive the exact trace formula for the Rydberg 
spectrum given in 

For the Hamiltonian we have derived a uniform trace formula for the diameter and circle orbit contributions 
that is valid for arbitrary e. We find that for low energies and low repetition numbers no tori exist, so that the gross- 
shell structure of this system is dominated by the circle and diameter orbits and their lowest repetitions. The tori 
bifurcate, at sufficiently high energy, out of the repetitions with k > 3 of the circle orbits, as was also observed recently 
for homogeneous power-law potentials in [32l |. Their contributions sufficiently high above the respective bifurcation 
energies are given by the usual Berry- Tabor type trace formula l|75|l . We have not discussed here the common uniform 
treatment of the circle orbits and the tori bifurcating from them, since this scenario is identical to the one discussed in 
detail in |3l| . The same qualitative results are found also in the limit e — * oo in which the potential becomes a purely 
quartic oscillator. Although the tori with N:M > 7:3 here exist at all energies, they are so much longer (by a factor 
> 10) than the shortest diameter and circle orbits that they only affect the finer details of the quantum spectrum. 

We thus have found the interesting and quite atypical situation of a three-dimensional system in which the periodic 
orbits with highest, i.e., three-fold degeneracy are only responsible for finer details of the quantum spectrum, whereas 
its gross-shell structure is to an astonishing degree dominated by the orbits of next-to- leading order in h, i.e., the 
diameter and circle orbits occurring in two-fold degenerate families. This is totally different from the situation, 
observed first by Balian and Bloch [J , of a spherical cavity potential (with ideally reflecting walls) which has turned 
out to be a realistic model for large spherical metal clusters. (We re-derive the semiclassical trace formula of Balian and 
Bloch for the spherical cavity in Appendix lD 11 ) There the famous super-shells (THESES come from the interference 
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of the shortest three-fold degenerate tori (the triangle and square orbits), whereas the two-fold degenerate diameter 
orbit family has virtually no effect on the observable shell structure of these systems. 

Our results explain qualitatively the numerical quantum-mechanical HF results [2lJ for harmonically trapped 
fermionic atom gases. A more quantitative comparison, in which the strength e of the perturbation in y| is de- 
termined directly from the numerically obtained self-consistent HF potentials, is in progress pjof. 

Lazzari et al. have studied the periodic orbits in a spherical Woods-Saxon potential in connection with the 
physics of metal clusters and thereby focused on the contributions of the diameter orbits. We find agreement of 
their results for the tori and the diameters with our results l|75|) . (|77|l and IjMfl . respectively, if we interpret their 
function t(Lm,E) as being half the radial period T r (E, Lnm) at the corresponding stationary value Lmi of the 
angular momentum (= for the diameters). They have, however, neglected the contributions from the circle orbits; 
furthermore they did not compare their semiclassical results with quantum-mechanical results. We suspect that the 
neglect of the circle orbits might be justifiable due to the chosen steepness of their Woods-Saxon potential, for which 
the triangle-like tori are perhaps more important than the circles. 

Ozorio de Almeida et al. j^l have discussed the summation of all periodic orbits in integrable systems in connection 
with spectral determinants. They conclude that edge corrections corresponding to lower-degenerate orbits may be 
neglected because they are implicitly approximated by the nearest-lying rational tori (characterised by very large num- 
bers M, N in our notation). While their argument applies to the full quantisation of the exact (or EBK- approximated) 
quantum spectrum obtained by summing over all orbits, it cannot be used when the spectrum is coarse-grained with a 
finite energy window 7, such as we have used the periodic-orbit sum in our present investigation (cf. also the discussion 
in the appendix of Ref. 351]) which is more practically oriented. 
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APPENDIX A: PARAMETRISATION THE MANIFOLD CP 2 

We reproduce here the parametrisation of CP 2 given in [2^]. We start with three complex numbers Z a 6 C: 

(Z 1 ,Z 2 ,Z 3 ) = { ni ,n2e iV2 ,n 3 e iV3 ), (Al) 

with n a € R, < n a < 1 (a = 1, 2, 3), and Vj € R, < Vj < 2tt (j = 2, 3), whereby the n a are restricted to 

n\ + n\ + n\ = 1 . (A2) 

The Z a define the complex projective space CP 2 . The n a are in one-to-one correspondence with the points on the 
first octant of the 2-sphere S 2 , and the Vj form a 2-torus. At the edges of the octant, the torus contracts to a circle 
or to a point, and one or both of the i/j are not defined. 

Physically, we interpret the real and negative imaginary parts of the Z a as coordinates and momenta, respectively, 
of a starting point (r ,p ) = ({r Q (0)}, {p a {0)}) 

r a (0) = dleZ a , p a (0) = -QmZ a , (A3) 

in the 6-dimensional phase space (r(i),p(t)) = ({_£«(£)}, {Pa(t)}) of the 3-dimensional isotropic harmonic oscillator. 
The Hamiltonian Q, which has U(3) symmetry |20l |. is given by 

ff(r >P ) = ±X> 2 r*+i£] =E. (A4) 

a=l 

The 27r-periodic solutions of the equations of motion are then given by 

r a (t) = ne(Z a e lult ), p a (t) = -Qm(Z Q e Mt ). (A5) 
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The equation (|A2|I defines the energy shell, so that the six quantities {u>r a ,p a } are points on a 5-sphere, S 5 (see also 
Appendix ITU below - ) . By taking the total energy to be E = 1/2 and setting u> = 1, we fix the radius of the 5-sphere 
to be unity. Taking out the phase factor e lujt takes us from S 5 to the projective space CP 2 . 

Using the Fubini-Study metric, the authors of [23| find the squared distance ds 2 between two points on CP 2 to be 

ds 2 = dn\ + dn 2 + dn 2 + n 2 (l — n\) dv\ + n 2 (l — n 2 ) dv 2 — 2n 2 n 2 d^dv^ . (A6) 

The first three terms correspond to the standard metric on S 2 , and the last three to the metric on a flat 2-torus whose 
shape depends on where we are on the octant. 

We now wish to re-parametrise the n a in terms of the two polar angles (■#, tp) on S 2 by writing 

n 1 =cos??, «2 — sin $ cos ip , n 3 = sin ?9 sin (p , (A7) 

with ip £ [0, 7r/2]. Collecting our four angles in a vector ip by defining 

ip = (i>i,ih,i/>3,il>i) = [®,<P,V2,va) , (A8) 
we can write the squared distance ds 2 as 

4 

ds 2 = 2J 9ij dipi dtpj , (A9) 

and find the elements gij of the metric tensor to be gn — 1, 1722 = sin 2 ??, 512 = <?2i = 0, 533 = sin 2 i9 cos 2 ^ (1 — 
sin 2 i9 cos 2 (p), 1744 = sin 2 ^sin 2 (/3 (1— sin 2 z9 sin 2 ^), and 334 = (743 = — sin 4, i9cos 2 ^sin 2 (/3. From this, we get the integration 
measure on CP 2 to be 

dQcp 2 = {det^gi^^^dipid^ ^"03 ^04 = sin 3, i9cos^(i^ cos ip sin <p dp dv^dv^ . (A-10) 
The integrated volume of CP 2 , consistent with a general formula for CP", is 

f tt 2 ir n 

tt CP 2 = / dfl CP 2 = — 4= Q C pn = — - . (All) 

/ 2 nl 



APPENDIX B: INTEGRATION OF THE MODULATION FACTOR OVER S' 

Here we give a strict proof of Eq. 121|) of the modulation factor M. for the anharmonically perturbed HO system 
The most straightforward paramctrisation of the six-dimensional phase space of the unperturbed HO is given 
in terms of the coordinate vector r(t) — {r a (t)} and the momentum vector p(t) = {p a (t)} (a = 1,2,3). Since a 
periodic solution is uniquely determined by their initial values at t — 0: r = r(0) and p = p(0), and the total energy 
E = (p 2 + lo 2 y 2 )/2 is a constant of motion, the six components of ro and po /u> must cover a 5-sphere with radius 
R = \/2E/ui. We introduce dimensionless coordinate and momentum vectors p and 7r, respectively, by 

P=^r , -=^Po, (Bl) 

so that in these variables the energy shell becomes the unit sphere S 5 : 

p 2 + TV 2 = 1 . (B2) 
From Eq. l|18fl . the first-order action perturbation ASi is given in terms of energy E and angular momentum L by 

ASx — — (7 (3 — lu 2 L 2 /E 2 ) , (B3) 
with a defined in (|19J) . The modulation factor for the first-order perturbed HO system becomes therefore 

M(ka/h) = ^ j dfl S 5 e lkAs ^ h , (B4) 



where f2g5 is the five-dimensional solid angle and dtts 5 is the integration measure on S 5 with J dfl s s = tt 3 . One 
may now express A5i in (|B3|I directly in terms of the five polar angles of six-dimensional hyperspherical coordinates, 
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which becomes a very complicated function, and integrate l(B4() over all five angles. We have done this numerically 
and verified that it yields exactly the same result as the numerical 4-dimensional CP 2 integration of 1(1211 with the 
phase function (fH^I . 

It is . ho wever, not necessary to perform the full five-dimensional integral for (|B4f) . We follow a much more economical 
route [33, exploiting the rotational symmetry of the system and the SO(3) invariance of the expression (|B3|I for A Si. 
Due to the restriction (|B2Jl . we can write the S 5 sphere as 

S 5 = j(costfe p ,sintf e7r ) € [0, tt/2] ; e p , e x e S 2 | , (B5) 

where e p and e,r are the unit vectors in the directions of p and 7r, respectively. The square of the conserved total 
angular momentum then is 
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L 2 = L 2 = fr x p ) 2 = —5- cos 2 tf sin 2 i? (e„ x ej 2 , (B6) 

so that the action perturbation becomes, after the substitution a = 2$, 

ASi = -er[3-sin 2 a(e p x e x ) 2 ] , ae[0,7r]. (B7) 
The integration measure for (|B5|) is 

dfl s s = cos 2 i9 dfl S 2(e p ) sm 2 i3 dfl S 2 (e n ) dd = - sin 2 a da dQ S 2 (e p ) dQ S 2 (e,) , (B8) 

8 

and the modulation factor becomes 

M(ka/h) = ^-3^ dasm 2 a J dO^e,) J dft S 2(e w ) e -*M3- 8 m 2 a (<V*^) 3 ]M (B9) 

We now introduce the angle 9 between the unit vectors e p and e T , so that 

(e p x ejr ) 2 = sin 2 6», (BIO) 

and choose e p as the direction of the north pole (9 = 0) of polar angles (9, <f>) for the vector e ff . The integrand then 
becomes independent of e p , so that the corresponding S 2 integral just gives J G?f2g2 (e p ) = An. For the other S 2 integral 
we have JdO,g2{e 7r ) — 2tt sin 9 d9 since the integrand does not depend on </>. We thus obtain, after the standard 
substitution u = cos 9 and an obvious reduction of the integration limits, 

M{ka/K) = - dasm 2 a due -^l3-sin 2 a (i-u 2 )]/h^ ( BU ) 
n Jo Jo 

Next, we make the substitution t = usina and then a further substitution s = cos a to obtain 

M(ka/h) = - da sina / dt e ^ka(3- si n 2 a +t 2 )/n , = f / ^ / ^ e _ ifcCT(2 + s ^)/fc (m2) 

t Jo Jo Jo Jo 

Since the last integral goes exactly over the first quadrant of a unit disk in the (s, t) plane, we can use polar coordinates 
(r, ip) and furthermore r 2 = 1 — z to obtain our final result 

M(ka/h) = - [ 1 dip ( drre- lk ^ 2+r ^/ h = / dz e -**«r(3-«)/Jl ) (B13) 



which is identical to that given in (|21() or in (|56(l . 



APPENDIX C: EXPLICIT INTEGRALS FOR THE ANHARMONICALLY PERTURBED HO 

In this appendix we derive analytical expressions for the particular scaled potential v(q) in 1)36(1 . Rewriting the 
radial action integral 1(43(1 using the substitution q 2 = z, we obtain 

P12 1 fZ2 J 

s r (e,l) = 2 d q ^/2e-q 2 -qy2-P/ q 2 = ^= / -V^), (CI) 
J 9 i V 2 J 2l 2 
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where 

g(z) = Aez - 2z 2 - z z - 2l 2 = (zi - z){z 2 - z)(z 3 - z) . (C2) 

The classical turning points are defined by the roots of the cubic equation g(z) = 0, which are always real: the 
discriminant can easily be shown to be negative except for I = Z m (e) where it becomes zero due to the double root 
z\ = z 2 . Selecting the roots such that 

z 3 < zi < z 2 , (C3) 

one finds that z 3 is always negative, whereas Z\ and z 2 are positive definite and represent the squares of the real 
classical turning points q\ and q%. 

The integrals in (|C1(I cannot be found in most tables. However, after an integration by parts we may write the 
r.h.s. as 

s r (e, I) = sP (e, I) + 4 2 ) (e, I) - fs^ (e, Z) , (C4) 

where we have defined 

S W( e ,Z) = V2/ dz^=. (C5) 
Jzi V9\ z ) 

In Byrd and Friedman |37| , these integrals are found to be: 

S «(e,Z) = 2 ^ ^[{n 2 -a 2 )K(n)+a 2 E{n)] 1 

4 2 )( e ,0 = 2 ^ A- [( 3w 4 _ 6a 2 K 2 + 2a 4 + k 2 q4) k{k) + 2a 2 ( g K 2 _ Q 2 _ ^2) E ^ 

V z 2 - z 3 3k 4 

S (-D( e ,0 = -^=ln( a 2 , K ). (C6) 

V z 2 - Z 3 Z 2 

Here K(k), E(k) and n(a 2 , k) are the complete elliptic integrals of first, second and third kind with modulus n, and 



Z2-Z1 2 1 Zx 



a 2 = 1 - =± . (C7) 

Z2 - Z 3 Z 2 

Note that < k 2 < a 2 < 1, which corresponds to the 'circular case' of the elliptic integral H(a 2 , k) |32|. This function 
becomes singular like (1 — a 2 ) -1 / 2 , i.e., for a — > 1 which happens when Z — > 0. Since zi goes to zero like l 2 /2e in this 
limit, n(a 2 , k) has a first-order pole at I = 0. Using the Laurent expansion of n(a 2 , k) given in Eq. 906.04 of [37|], we 
find 

(e, I) = - ^ - 2a(e) + 0(f) = ~ - 2o(e) + O(0 , (C8) 

y/~Z 1 Z 2 Z 3 I 



with 



a{e) = 



. E(k ) 



i('-T' + -l- 



Here z z ° and «o are the f = values of the quantities given in Eqs. (|CJ14|) and l|C15fl below. The last equality in (|C8|) 
follows because z±z 2 z 3 — —2l 2 , as seen from (IC2|) . The functions Sj- (e, Z) with n = 1,2 are both regular in Z = 0. 
Altogether we obtain the following Taylor expansion for the total radial action integral l|C4() : 



s r (e, = s r (e, 0) - 7T Z + <z(e) Z 2 + C(Z 3 ) , (C10) 
The radial period (|44|l is easily found to be 

f r (e,Z) = 4°)(e,Z)= 2v/ ^ (Cll) 

V Z 2 — ^3 
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In the limit I = 0, which is relevant for the diameter orbits and the TF density of states, the above results simplify 
considerably. We obtain 



*r(e,0) = - (1 + Ae) 1 /* [(V (1 + 4e) + l) K(k ) - 2E(k )] , (C12) 
tr(e,0) = 2{l + 4e)- 1/4 K(n ) . (C13) 

Note that in this limit we get 

z° = G, z° = 0?°) 2 = VT+4e- 1, ^ = -VI + 4e - 1 , (C14) 

and the modulus of the elliptic integrals becomes 



2VTT4^ 

The TF expression (|48|) for the density of states becomes 

STF(e) = — ^— ^J(l + 4e) 1 / 4 [2(l + 3e)£;( Ko )-(l + 3e + VTT4^)X( K o)]. (C16) 

It is of interest to see the first terms in the Taylor expansion of the above results around e = 0, whose leading terms 
are those of the pure harmonic oscillator: 

s r (e, /) = 7T (e - - - 7T (3e 2 - I 2 ) + . . . , (C17) 
8 

t r (e,l) = 7r-^7re+^7r(7e 2 -Z 2 ) + ... , (C18) 
4 64 

5TF(e) = ^L^e 2 (l-5e/4 +...). (C19) 

The unsealed leading HO terms are 

S?(E,L)=J^-l\=it[LM-L] ) 7?° = -, gfeffl = -&. (C20) 

The diameter orbit with L = makes two full radial oscillations during its primitive period. Therefore its primitive 
action and period are given by 

S d (e) = 2S r (e, 0) , T d {e) = 2T r (e, 0) . (C21) 

For the circle orbit, the primitive action integral S c (e) is just 2tt times the maximum value of the angular momentum, 
L ra (e) = sl m (e), which for the present potential is found to be 



L m (e) = sJ J- [(1 + 3e)3/ 2 - 1 - 9e/2] , S c (e) = 2vL m {e) = — (1 - e/4 + . . . ) 



(C22) 



The primitive period of the circle orbit is given by its energy derivative 



T c (e) = ^ - (l + ^-l = *L (1 _ e/2 + . . . ) . (C23) 
w v/(l + 3e) 3 / 2 - 1 - 9e/2 w 

We also give here the value of the radial period T r (e,l) at Z = i m (e). The corresponding frequency oj r (e,l m ) = 
2ir /T r (e,l m ) is easily obtained from the second derivative of the effective potential at its minimum qoj which for the 
present potential is given by 

g 2 = z = 2(VTT3i-l)/3, (C24) 
using uj 2 (e,l m ) — V^'fAqo), and the result becomes 

T r (e, U = — -r— ; — r = - (1 + ie)- 1 ^ = - (1 - 3e/4 +...). (C25) 
Lu r (e,l m ) oj uo 

This agrees with the result obtained from i|Cll|) . noting that for the maximum value of I = l m (e), Z\ = z 2 = z and 
z 3 = -2(2yT + 3e + l)/3, so that k JCTJ becomes zero and K(0) = tt/2. 
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APPENDIX D: TRACE FORMULAE FOR OTHER SPHERICAL POTENTIALS 

In order to illustrate the validity and the usage of our general EBK trace integral (|46|l . we present in this appendix 
briefly its application to two popular spherical potentials. 

1. The spherical box potential 

We consider a particle with mass m = 1 in a spherical box with radius R and ideal specular reflection from the 
boundary i.e., a three-dimensional spherical billiard. Quantum- mechanically, one obtains the spectrum by solving the 
Schrodinger equation with Dirichlet boundary conditions. The semiclassical trace formula for this system has been 
derived by Balian and Bloch Q using a multi-reflection expansion of the Green function. To derive it from our formula 
H46JI . we need the radial action integral for arbitrary energy E > and angular momentum < L < L m — R\j2E. 
Working throughout with unsealed variables, we find (with r% = L/\/2E) 

rR 

S r (E,L) = 2 / dry/2E - L 2 /r 2 = 2y / 2ER 2 - L 2 -2L arccos (L/RV2E) 

= 2RV2E-t:L + L 2 /rV2E + 0(L 3 ), (Dl) 

thus fulfilling our general relation (|C10|) with 

a(e) = 1/RV2E. (D2) 

With this, the stationary condition (|54|l gives 

2N arccos (L/RV2E) = 2nM (D3) 

with the solutions 

L NM = RV2Ecos(ttM/N). (D4) 



As in the potentials discussed in Sect. IIIII the solution with L = corresponding to N=2M contributes only via the 
end-point correction of (|46|l to the diameter orbit with the primitive action 

S d (E) = 2S r {E, 0) = ARV2E (D5) 
and, according to the lower equation in (|72|l . with the amplitude of its fc-th repetition 

Ai(E) = ~. (k = \M\) (D6) 

The upper end-point correction with Lnm = L m = R\f2E can, for M ^ 0, only be reached formally in the limit 
N — > oo. It corresponds to the "whispering gallery mode" with amplitude A oc 1/yN, cf. (|D8|I below, and does 
therefore not contribute to the trace formula. The tori with < Lnm < L m have the actions 

Sl M {E) = NS r (E, L NM ) + 2nML NM = 2NRV2E sin <p NM , ?nm = kM/N . (D7) 



' NM 

This corresponds exactly to the actions of the polygonal orbits given in with winding number fc = \M\ > 1 and 
v = \N\ > 2k corners (i.e., reflections from the boundary). The case v = 2k corresponds to the k-th repetition of 
the diameter orbit. For the polygons, cp v k = nk/v is half the polar angle covered by one of their segments. The 
amplitudes of the tori with v > 2k become, according to (JHJ, 



, T , . 2R b / 2 (2E) 1 / A , . /sin w vk , . 

A ^ E ) - hB /2j M*Pv*))l ■ (D8) 

Before we can write down the trace formula with the correct phases, we have to correct the radial quantisation 
condition (|27|) . because the Maslov index changes from 1/4 to 1/2 per turning point if a reflection happens there. 
Hence the quantisation condition is S r = 27th(n r + 1) for the diameter orbit and S r = 2irh(n r + 3/4) for the tori. This 
changes the phases in (|46H and all the results derived from it. Taking this into account, we obtain the semiclassical 
trace formula for the spherical billiard 

Sg(E) = J2 A vk(E) sm[S^ k {E)/h-3vTr/2~kn + 3tt/4] - Y, A k( E ) sin[kS d {E)/h] , (D9) 

v>2k>0 k>Q 

which is identically the result of Balian and Bloch 2] . 
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2. The Coulomb potential 

Here we derive the trace formula for the Coulomb potential as a special central potential which is not regular in 



0: 



V(r) = - 



(D10) 



We need not use scaled variables here, but it is useful to introduce the positive energy e = — E since all bound states 
in the potential (|D10fl have E < 0. The maximum angular momentum and the action of the circle orbit are 



L m (e) 



2e 



S c (e) 



2na 



2e 



from which the period of the circle orbit is found to be 

_dL m (e) 



T c {e) = 2vr- 



dE s/2e 3 / 2 
The radial action integral becomes elementary (cf. 

a 



rr 2 

S r (e, L) = V2 dr^J-e + a/r - L 2 /2r 2 
Jn 

From it we find the radial period 



2tt 



2e 



L ) =27r[i ro (e)-i] 



T r (e) 



V2e 3 / 2 



(Dll) 



(D12) 



(D13) 



(D14) 



The difference from the other potentials treated in this paper is that the term linear in L in (|D13(1 here has the 
coefficient — 2ir. Furthermore, we find from the general relation 1)521) 



dS r (e,L) 



= -2tt 



u^(e) = w r (e) , T r (e) = T^(e) . 



(D15) 



Thus, the radial and angular frequencies arc identically the same for all values of the angular momentum L. This is 
a special property of Kepler's ellipses. [Note that here, in contrast to the spherical HO, the angular momentum L is 
defined with respect to one of the focal points and not to the symmetry centre of the ellipse!] 

Inserting the above results into the general EBK trace integral l|46[) for the density of states, we obtain 



fa(e) = f £ (-l) M+N 



irh 3 



M,N=-oo 



LdL e i2w[N(L m -L)+ML]/H_ 



(D16) 



Due to (|D15|) . the general resonance condition (|54|) implies N = M for all periodic orbits. We therefore only use the 
terms with N = M = ±k in the double sum of l|D16|) (cf. the discussion at the end of Sect. IIII . We then find 



9coui{e) 



h 3 (2e) 5 / 2 



1 + 2J2 



k 27ra\ 
cos I — 



fe=l 



Writing this in terms of atomic units (note that we have put m = 1) 

Eat = 2iP =Ry > 

we finally obtain the quantum-mechanically exact trace formula 



Qcoul 



{E) = Y^n 2 8{E + E at /n 2 ) = 



e; 



3/2 



2(-Efn 



oo 

1 + 2 E cos (2irkyJ-E at /E 



k=l 



which has been given in The eigenenergies 



E n = 



E a 



(n = n r + £+ 1) 



(D17) 



(D18) 



(D19) 



(D20) 



with degeneracies d n = n 2 are, as is well known, also found directly from the EBK quantisation conditions l|27(l - l|29() 
using |D13)l and 
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